The likelihood function

The maximum likelihood estimate

Finding the maximum of a function

Recall from calculus: the maxima of a continuous function \(f(x)\) over a range of input values \(\mathcal{R}(x)\) can be found either:

• at the edges of the range of input values, OR:
• where the function is flat (i.e. where the gradient function \(f'(x) = 0\)) AND the second derivative is negative definite (\(f''(x) < 0\)).

Directly maximizing the likelihood function for independent data

To find the maximizer(s) of the likelihood function, we need to solve \(\mathscr{L}'(\theta) = 0\) for \(\theta\). However, even for mutually independent data, we quickly run into a problem:

\[ \begin{aligned} \mathscr{L}'(\theta) &= \frac{\partial}{\partial \theta} \mathscr{L}(\theta) \\ &= \frac{\partial}{\partial \theta} \prod_{i=1}^n p(X_i=x_i|\theta) \end{aligned} \tag{3}\]

The derivative of the likelihood of independent data is the derivative of a product. To evaluate this derivative, we will have to perform a massive application of the product rule (?@thm-product-rule).

\[\ell'= \frac{\varepsilon}{\text{Var}{\left(X\right)}}\]

Asymptotic distribution of the maximum likelihood estimate

Recall:

We can estimate \(\mathcal{I}(\theta)\) using either \(\mathcal{I}(\hat\theta_{\text{ML}})\) or \(I(\tilde{x}; \hat\theta_{\text{ML}})\).

So we can estimate the standard error of \(\hat\theta_k\) as:

\[ \widehat{\text{SE}}{\left(\hat\theta_k\right)} = \sqrt{{\left[{\left(\hat{\mathcal{I}}{\left(\hat\theta_{\text{ML}}\right)}\right)}^{-1}\right]}_{kk}} \]

The (Fisher) (expected) information matrix

The variance of \(\ell'(x,\theta)\), \({Cov}\left\{ \ell'(x,\theta) \right\}\), is also very important; we call it the “expected information matrix”, “Fisher information matrix”, or just “information matrix”, and we represent it using the symbol \(\mathcal{I}{\left(I\right)}\) (\scriptI in Unicode, \mathcal{I} in LaTeX).

\[ \begin{aligned} \mathcal{I} \stackrel{\text{def}}{=}\mathcal{I}(\theta) \\ &\stackrel{\text{def}}{=}\text{Cov}{\left(\ell'|\theta\right)} \\ &= \text{E}[ \ell'{\ell'}^{\top} ] - \text{E}[\ell'] \ \text{E}[\ell']^{\top} \end{aligned} \]

The elements of \(\mathcal{I}\) are:

\[ \begin{aligned} \mathcal{I}_{ij} &\stackrel{\text{def}}{=}\text{Cov}{\left({\ell'}_{i},{\ell'}_{j}\right)} \\ &= \text{E}[ \ell_{i}'\ell_{j}' ] - \text{E}[ {\ell'}_{i} ] \text{E}[ {\ell'}_{j} ] \end{aligned} \]

Here,

\[ \begin{aligned} \text{E}{\left[\ell'\right]} &\stackrel{\text{def}}{=}\int_{x \in \mathcal{R}(x)} { \ell'(x,\theta) \text{p}(X = x | \theta) dx } \\ &= \int_{x \in \mathcal{R}(X)} { {\left( \frac{\partial}{\partial \theta} \log{\text{p}(X = x | \theta)} \right)} \text{p}(X = x | \theta) dx } \\ &= \int_{x \in \mathcal{R}(X)} { \frac {\frac{\partial}{\partial \theta} \text{p}(X = x | \theta)} {\text{p}(X = x | \theta)} \text{p}(X = x | \theta) dx } \\ &= \int_{x \in \mathcal{R}(X)} { \frac{\partial}{\partial \theta} \text{p}(X = x | \theta) dx } \end{aligned} \]

And similarly

\[ \text{E}{\left[\ell' \ell'^{\top}\right]} \stackrel{\text{def}}{=} \int_{x \in R(x)} {\ell'(x,\theta)\ell'(x,\theta)^{\top}\ \text{p}{\left(X = x | \theta\right)}\ dx} \]

Note that \(\text{E}{\left[\ell'\right]}\) and \(\text{E}{\left[\ell'{\ell'}^{\top}\right]}\) are functions of \(\theta\) but not of \(x\); the expectation operator removed \(x\).

Also note that for most of the distributions you are familiar with (including Gaussian, binomial, Poisson, exponential):

\[\text{E}{\left[\ell'\right]} = 0\]

So

\[\mathcal{I}{\left(\theta\right)} = \text{E}{\left[\ell'{\ell'}^{\top} \right]}\]

Moreover, for those distributions (called the “exponential family”), we have:

\[ \mathcal{I}= -\text{E}{\left[\ell''\right]} = \text{E}{\left[- \ell'\ell'^{\top}\right]} \]

(see Dobson and Barnett (2018), §3.17).

Observed information

Sometimes, we use \(I(\theta;x) \stackrel{\text{def}}{=}- hess\) (note the standard-font “I” here). \(I(\theta;x)\) is the observed information, precision, or concentration matrix (Negative Hessian).

We can estimate \(\mathcal{I}^{- 1}(\theta)\) by working out \(\text{E}{\left[-\ell''\right]}\) or \(\text{E}{\left[\ell'{\ell'}^{\top}\right]}\) and plugging in \(\hat\theta_{\text{ML}}\), but sometimes we instead use \(I(\hat\theta_{\text{ML}},\tilde{x})\) for convenience; there are some cases where it’s provably better according to some criteria (Efron and Hinkley (1978)).

Quantifying (un)certainty of MLEs

Confidence intervals for MLEs

An asymptotic approximation of a 95% confidence interval for \(\theta_k\) is

\[ \hat\theta_{\text{ML}}\pm z_{0.975} \times \widehat{\text{SE}}{\left(\hat\theta_k\right)} \]

where \(z_\beta\) the \(\beta\) quantile of the standard Gaussian distribution.

p-values and hypothesis tests for MLEs

(to add)

Likelihood ratio tests for MLEs

log(likelihood ratio) tests (c.f. Dobson and Barnett 2018, sec. 5.7):

\[2{\left(\ell- \ell_{\text{H0}}\right)} \sim \chi^{2}(p - q)\]

See also https://online.stat.psu.edu/stat504/book/export/html/657

Prediction intervals for MLEs

\[\overline{X} \in [ \widehat{\mu} \pm z_{1 - \alpha\text{/}2}\frac{\sigma}{m} ]\]

Where \(m\) is the sample size of the new data to be predicted (typically 1, except for binary outcomes, where it needs to be bigger for prediction intervals to make sense).

Data

The cyclones dataset in the dobson package (Table 1) records the number of tropical cyclones in Northeastern Australia during 13 November-to-April cyclone seasons (more details in Dobson and Barnett (2018) §1.6.5 and help(cyclones, package = "dobson")). Figure 1 graphs the number of cyclones (y-axis) by season (x-axis). Let’s use \(Y_i\) to represent these counts, where \(i\) is an indexing variable for the seasons and \(Y_i\) is the number of cyclones in season \(i\).

Exploratory analysis

Suppose we want to learn about how many cyclones to expect per season.

library(dobson)
library(dplyr)
data(cyclones)
library(pander)
pander(cyclones |> relocate(season, .before = everything()))

library(ggplot2)
library(dplyr)
cyclones |>
  mutate(years = factor(years, levels = years)) |>
  ggplot(aes(x = years, y = number, group = 1)) +
  geom_point() +
  geom_line() +
  xlab("Season") +
  ylab("Number of cyclones") +
  expand_limits(y = 0) +
  theme(axis.text.x = element_text(vjust = .5, angle = 45))

There’s no obvious correlation between adjacent seasons, so let’s assume that each season is independent of the others.

Let’s also assume that they are identically distributed; let’s denote this distribution as \(P(Y=y)\). Note that there’s no index \(i\) in this expression, since we are assuming the \(Y_i\)s are identically distributed.

We can visualize the distribution using a bar plot (Figure 2).

[R code]

cyclones |>
  ggplot() +
  geom_histogram(aes(x = number)) +
  expand_limits(x = 0) +
  xlab("Number of cyclones") +
  ylab("Count (number of seasons)")

Table 2 provides summary statistics.

n <- nrow(cyclones)
sumx <- cyclones |>
  pull(number) |>
  sum()
xbar <- cyclones |>
  pull(number) |>
  mean()

cyclones |> table1::table1(x = ~number)

Model

We want to estimate \(P(Y=y)\); that is, \(P(Y=y)\) is our estimand.

We could estimate \(P(Y=y)\) for each value of \(y\) in \(0:\infty\) separately (“nonparametrically”) using the fraction of our data with \(Y_i=y\), but then we would be estimating an infinitely large set of parameters, and we would have low precision. We will probably do better with a parametric model.

Estimating the model parameters using maximum likelihood

Now, we can estimate the parameter \(\lambda\) for this distribution using maximum likelihood estimation.

The score function

The Hessian matrix

Finding the MLE analytically

In this case, we can find the MLE of \(\lambda\) by solving the score equation for \(\lambda\) analytically (using algebra):

Let’s call this solution of the estimating equation \(\tilde \lambda\) for now:

\[\tilde \lambda \stackrel{\text{def}}{=}\bar x\]

Information matrices

[R code]

obs_inf <- function(...) -hessian(...)
ggplot() +
  geom_function(fun = obs_inf, n = 1001) +
  xlim(min(cyclones$number), max(cyclones$number)) +
  ylab("I(lambda)") +
  xlab("lambda") +
  geom_hline(yintercept = 0, col = "red")

Iterative maximization

(c.f., Dobson and Barnett (2018), Chapter 4)

\[ \begin{aligned} {\widehat{\theta}}^* &\leftarrow {\widehat{\theta}}^* + {\left(I{\left(\tilde{y};{\widehat{\theta}}^*\right)}\right)}^{-1} \ell'{\left(\tilde{y};{\widehat{\theta}}^*\right)} \\ &= {\widehat{\theta}}^* - {\left(\ell''{\left(\tilde{y};{\widehat{\theta}}^*\right)}\right)}^{-1} \ell'{\left(\tilde{y};{\widehat{\theta}}^*\right)} \end{aligned} \]

\[ \begin{aligned} \ell'(\theta) &\approx \ell'^*(\theta) \\ &\stackrel{\text{def}}{=}\ell'({\widehat{\theta}}^*) + \ell''({\widehat{\theta}}^*)(\theta- {\widehat{\theta}}^*) \end{aligned} \]

\[ \begin{aligned} \theta &= {\widehat{\theta}}^* - \ell'({\widehat{\theta}}^*) \cdot{\left(\ell''({\widehat{\theta}}^*)\right)}^{-1} \end{aligned} \]

For computational simplicity, we will sometimes use \(\mathfrak{I}^{- 1}(\theta)\) in place of \(I\left( \widehat{\theta},y \right)\); doing so is called “Fisher scoring” or the “method of scoring”. Note: this substitution is the opposite of the substitution that we are making for estimating the variance of the MLE; this time we should technically use the observed information but we use the expected information instead.

There’s also an “empirical information matrix” (see McLachlan and Krishnan (2007)):

\[I_{e}(\theta,y) \stackrel{\text{def}}{=}\sum_{i = 1}^{n}{\ell_{i}'\ {\ell_{i}'}^{\top}} - \frac{1}{n}\ell'{\ell'}^{\top}\]

where \(\ell_{i}\) is the log-likelihood of the ith observation. Note that \(\ell' = \sum_{i = 1}^{n}\ell_{i}'\).

\(\frac{1}{n}I_{e}(\theta,y)\) is the sample equivalent of

\[\mathfrak{I \stackrel{\text{def}}{=}I(}\theta) \stackrel{\text{def}}{=}{Cov}\left( \ell'|\theta \right) = E[ \ell'{\ell'}^{\top} ] - E[ \ell' ]\ E[ \ell' ]^{\top}\]

\[\left\{ \mathfrak{I}_{jk} \stackrel{\text{def}}{=}{Cov}\left( {\ell'}_{j},{\ell'}_{k} \right) = E[ \ell_{j}'\ell_{k}' ] - E[ {\ell'}_{j} ] E[ {\ell'}_{k} ] \right\}\]

\(I_{e}(\theta,y)\) is sometimes computationally easier to compute for Newton-Raphson-type maximization algorithms.

c.f. https://en.wikipedia.org/wiki/Newton%27s_method_in_optimization

\[ \ell'( \lambda; \tilde{x} ) = {\left(\frac{72}{\lambda} \right)} - n \]

\[ \ell''( \lambda; \tilde{x} ) = -\frac{72}{\lambda^2} \]

\[ \begin{aligned} \ell'(\lambda) &\approx \ell'^*(\lambda) \\ &\stackrel{\text{def}}{=}\ell'({\widehat{\lambda}}^*) + \ell''({\widehat{\lambda}}^*)(\lambda - {\widehat{\lambda}}^*) \\ &= {\left(\frac{72}{{\widehat{\lambda}}^*} - n\right)} + {\left(-\frac{72}{{\left({\widehat{\lambda}}^*\right)}^2}\right)} (\lambda - {\widehat{\lambda}}^*) \end{aligned} \]

\[ \ell^*(\lambda) = \ell({\widehat{\lambda}}^*) + (\lambda - {\widehat{\lambda}}^*) \ell'({\widehat{\lambda}}^*) + \frac{1}{2}\ell''({\widehat{\lambda}}^*)(\lambda-{\widehat{\lambda}}^*)^2 \]

approx_loglik <- function(lambda, lhat, ...) {
  loglik(lambda = lhat, ...) +
    score(lambda = lhat, ...) * (lambda - lhat) +
    1 / 2 * hessian(lambda = lhat, ...) * (lambda - lhat)^2
}

plot_loglik <- ggplot() +
  geom_function(
    fun = loglik,
    aes(col = "true log-likelihood"),
    n = 1001
  ) +
  geom_function(
    fun = approx_loglik,
    aes(col = "approximate log-likelihood"),
    n = 1001,
    args = list(lhat = cur_lambda_est)
  ) +
  geom_point(
    size = point_size,
    aes(
      x = cur_lambda_est, y = loglik(lambda = cur_lambda_est),
      col = "current estimate"
    )
  ) +
  geom_point(
    size = point_size,
    aes(
      x = xbar,
      y = loglik(xbar),
      col = "true MLE"
    )
  ) +
  xlim(min(cyclones$number) - 1, max(cyclones$number)) +
  ylab("l'(lambda)") +
  xlab("lambda")

print(plot_loglik)

\[ \begin{aligned} \lambda &= {\widehat{\lambda}}^*- \ell'({\widehat{\lambda}}^*) \cdot{\left(\ell''({\widehat{\lambda}}^*)\right)}^{-1} \\ &= 4.375 \end{aligned} \]

new_lambda_est <-
  cur_lambda_est -
  score(cur_lambda_est) * hessian(cur_lambda_est)^-1

The score function

\[\ell'(x,\theta) \stackrel{\text{def}}{=}\frac{\partial}{\partial \theta}\ell(x,\theta) = \left( \begin{array}{r} \frac{\partial}{\partial \mu}\ell(\theta;x) \\ \frac{\partial}{\partial \sigma^{2}}\ell(\theta;x) \end{array} \right) = \left( \begin{array}{r} \ell_{\mu}'(\theta;x) \\ \ell_{\sigma^{2}}'(\theta;x) \end{array} \right)\].

\(\ell'(x,\theta)\) is the function we set equal to 0 and solve to find the MLE:

\[{\widehat{\theta}}_{ML} = \left\{ \theta:\ell'(x,\theta) = 0 \right\}\]

MLE of \(\mu\)

\[ \begin{aligned} \frac{d\ell}{d\mu} &= - \frac{1}{2}\sum_{i = 1}^{n} \frac{- 2(x_{i} - \mu)}{\sigma^{2}} \\ &= \frac{1}{\sigma^{2}} {\left[ {\left( \sum_{i = 1}^{n}x_{i} \right)} - n\mu \right]} \end{aligned} \]

If \(\frac{d\ell}{d\mu} = 0\), then \(\mu = \overline{x} \stackrel{\text{def}}{=}\frac{1}{n}\sum_{i = 1}^{n}x_{i}\).

\[\frac{d^{2}\ell}{(d\mu)^{2}} = \frac{- n}{\sigma^{2}} < 0\]

So \({\widehat{\mu}}_{ML} = \overline{x}\).

MLE of \(\sigma^{2}\)

\[\frac{d\ell}{d\sigma^{2}} = \frac{\partial}{\partial \sigma^{2}}\left( - \frac{n}{2}\log{\sigma^{2}} - \frac{1}{2}\sum_{i = 1}^{n}\frac{\left( x_{i} - \mu \right)^{2}}{\sigma^{2}} \right)\ \]

\[= - \frac{n}{2}\left( \sigma^{2} \right)^{- 1} + \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2}\]

If \(\frac{d\ell}{d\sigma^{2}} = 0\), then:

\[\frac{n}{2}\left( \sigma^{2} \right)^{- 1} = \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2}\]

\[\sigma^{2} = \frac{1}{n}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2}\]

We plug in \({\widehat{\mu}}_{ML} = \overline{x}\) to maximize globally (a technique called profiling):

\[\hat{\sigma}_{ML}^{2} = \frac{1}{n}\sum_{i = 1}^{n}\left( x_{i} - \overline{x} \right)^{2}\]

Now:

\[ \begin{aligned} \frac{d^{2}\ell}{\left( d\sigma^{2} \right)^{2}} &= \frac{\partial}{\partial \sigma^{2}}\left\{ - \frac{n}{2}\left( \sigma^{2} \right)^{- 1} + \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2} \right\}\\ &= \left\{ - \frac{n}{2}\frac{\partial}{\partial \sigma^{2}}\left( \sigma^{2} \right)^{- 1} + \frac{1}{2}\frac{\partial}{\partial \sigma^{2}}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2} \right\}\\ &= \left\{ \frac{n}{2}\left( \sigma^{2} \right)^{- 2} - \left( \sigma^{2} \right)^{- 3}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2} \right\}\\ &= \left( \sigma^{2} \right)^{- 2}\left\{ \frac{n}{2} - \left( \sigma^{2} \right)^{- 1}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2} \right\} \end{aligned} \]

Evaluated at \(\mu = \overline{x},\sigma^{2} = \frac{1}{n}\sum_{i = 1}^{n}\left( x_{i} - \overline{x} \right)^{2}\), we have:

\[ \begin{aligned} \frac{d^{2}\ell}{\left( d\sigma^{2} \right)^{2}} &= \left( {\widehat{\sigma}}^{2} \right)^{- 2}\left\{ \frac{n}{2} - \left( {\widehat{\sigma}}^{2} \right)^{- 1}\sum_{i = 1}^{n}\left( x_{i} - \overline{x} \right)^{2} \right\}\\ &= \left( {\widehat{\sigma}}^{2} \right)^{- 2}\left\{ \frac{n}{2} - \left( {\widehat{\sigma}}^{2} \right)^{- 1}n{\widehat{\sigma}}^{2} \right\}\\ &= \left( {\widehat{\sigma}}^{2} \right)^{- 2}\left\{ \frac{n}{2} - n \right\}\\ &= \left( {\widehat{\sigma}}^{2} \right)^{- 2}n\left\{ \frac{1}{2} - 1 \right\}\\ &= \left( {\widehat{\sigma}}^{2} \right)^{- 2}n\left( - \frac{1}{2} \right) < 0 \end{aligned} \]

Finally, we have:

\[ \begin{aligned} \frac{d^{2}\ell}{d\mu\ d\sigma^{2}} &= \frac{\partial}{\partial \mu}\left\{ - \frac{n}{2}\left( \sigma^{2} \right)^{- 1} + \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2} \right\}\\ &= \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\frac{\partial}{\partial \mu}\sum_{i = 1}^{n}\left( x_{i} - \mu \right)^{2}\\ &= \frac{1}{2}\left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}{- 2(x_{i} - \mu)}\\ &= - \left( \sigma^{2} \right)^{- 2}\sum_{i = 1}^{n}{(x_{i} - \mu)} \end{aligned} \]

Evaluated at \(\mu = \widehat{\mu} = \overline{x},\sigma^{2} = {\widehat{\sigma}}^{2} = \frac{1}{n}\sum_{i = 1}^{n}\left( x_{i} - \overline{x} \right)^{2}\), we have:

\[\frac{d^{2}\ell}{d\mu\ d\sigma^{2}} = - \left( {\widehat{\sigma}}^{2} \right)^{- 2}\left( n\overline{x} - n\overline{x} \right) = 0\]

Covariance matrix

\[I = \begin{bmatrix} \frac{n}{\widehat{\sigma}^{2}} & 0 \\ 0 & \left( {\widehat{\sigma}}^{2} \right)^{- 2}n\left(\frac{1}{2} \right) \end{bmatrix} = \begin{bmatrix} a & 0 \\ 0 & d \end{bmatrix}\]

So:

\[I^{- 1} = \frac{1}{ad}\begin{bmatrix} d & 0 \\ 0 & a \end{bmatrix} = \begin{bmatrix} \frac{1}{a} & 0 \\ 0 & \frac{1}{d} \end{bmatrix}\]

\[I^{- 1} = \begin{bmatrix} \frac{{\widehat{\sigma}}^{2}}{n} & 0 \\ 0 & \frac{{2\left( {\widehat{\sigma}}^{2} \right)}^{2}}{n} \end{bmatrix}\]

See Casella and Berger (2002) p322, example 7.2.12.

To prove it’s a maximum, we need:

• \(\ell' = 0\)

• At least one diagonal element of \(\ell''\) is negative.

• Determinant of \(\ell''\) is positive.

n_obs <- 100 # we're going to take a small subset of the data to look at;
# if we took the whole data set, the likelihood function would be hard to
# graph nicely

library(dplyr)
data1 <-
  hers |>
  filter(
    diabetes == 0,
    exercise == 0
  ) |>
  head(n_obs)

glucose_data <-
  data1 |>
  pull(glucose)

library(ggplot2)
plot1 <-
  data1 |>
  ggplot(aes(x = glucose)) +
  geom_histogram(aes(x = glucose, after_stat(density))) +
  theme_classic()

print(plot1)

Looks somewhat plausibly Gaussian. Good enough for this example!

Find the MLEs

mu_hat <- mean(glucose_data)
sigma_sq_hat <- mean((glucose_data - mean(glucose_data))^2)

Our MLEs are:

• \(\hat\mu = 98.66\)

• \(\hat\sigma^2 = 104.7444\)

Here’s the estimated distribution, superimposed on our histogram:

plot1 +
  geom_function(
    fun = function(x) dnorm(x, mean = mu_hat, sd = sqrt(sigma_sq_hat)),
    col = "red"
  )

Looks like a somewhat decent fit? We could probably do better, but that’s for another time.

Construct the likelihood and log-likelihood functions

loglik <- function(
    mu, # I'm assigning default values, which the function will use
    # unless we tell it otherwise
    sigma = sd(x), # note that you can define some default inputs
    # based on other arguments
    x = glucose_data,
    n = length(x)) {
  normalizing_constants <- -n / 2 * log((sigma^2) * 2 * pi)

  likelihood_kernel <- -1 / (2 * sigma^2) * {
    # I have to do this part in a somewhat complicated way
    # so that we can pass in vectors of possible values of mu
    # and get the likelihood for each value;
    # for the binomial case it's easier
    sum(x^2) - 2 * sum(x) * mu + n * mu^2
  }

  answer <- normalizing_constants + likelihood_kernel

  return(answer)
}

# `...` means pass any inputs to lik() along to loglik()
lik <- function(...) exp(loglik(...))

Graph the Likelihood as a function of \(\mu\)

(fixing \(\sigma^2\) at \(\hat\sigma^2 = 104.7444\))

ggplot() +
  geom_function(fun = function(x) lik(mu = x, sigma = sigma_sq_hat)) +
  xlim(mean(glucose_data) + c(-1, 1) * sd(glucose_data)) +
  xlab("possible values of mu") +
  ylab("likelihood") +
  geom_vline(xintercept = mean(glucose_data), col = "red")

Graph the Log-likelihood as a function of \(\mu\)

(fixing \(\sigma^2\) at \(\hat\sigma^2 = 104.7444\))

ggplot() +
  geom_function(fun = function(x) loglik(mu = x, sigma = sigma_sq_hat)) +
  xlim(mean(glucose_data) + c(-1, 1) * sd(glucose_data)) +
  xlab("possible values of mu") +
  ylab("log(likelihood)") +
  geom_vline(xintercept = mean(glucose_data), col = "red")

Likelihood and log-likelihood for \(\sigma\), conditional on \(\mu = \hat\mu\):

ggplot() +
  geom_function(fun = function(x) lik(sigma = x, mu = mean(glucose_data))) +
  xlim(sd(glucose_data) * c(.9, 1.1)) +
  geom_vline(
    xintercept = sd(glucose_data) * sqrt(n_obs - 1) / sqrt(n_obs),
    col = "red"
  ) +
  xlab("possible values for sigma") +
  ylab("Likelihood")

ggplot() +
  geom_function(
    fun = function(x) loglik(sigma = x, mu = mean(glucose_data))
  ) +
  xlim(sd(glucose_data) * c(0.9, 1.1)) +
  geom_vline(
    xintercept =
      sd(glucose_data) * sqrt(n_obs - 1) / sqrt(n_obs),
    col = "red"
  ) +
  xlab("possible values for sigma") +
  ylab("log(likelihood)")

Standard errors by sample size:

Recall from Section 1.4.4 that the asymptotic standard error of \({\widehat{\mu}}_{ML}\) is

\[ \begin{aligned} \widehat{\text{SE}}{\left({\widehat{\mu}}\right)} &= \sqrt{{\left[{\left(\hat{\mathcal{I}}{\left(\widehat{\mu}_{ML}\right)}\right)}^{-1}\right]}} \\ &= \frac{\hat{\sigma}}{\sqrt{n}} \end{aligned} \]

se_mu_hat <- function(n, sigma = sd(glucose_data)) sigma / sqrt(n)
ggplot() +
  geom_function(fun = se_mu_hat) +
  scale_x_log10(
    limits = c(10, 10^5), name = "Sample size",
    labels = scales::label_comma()
  ) +
  ylab("Standard error of mu (mg/dl)") +
  theme_classic()

Power

Rejection region

For example, suppose we wish to detect a difference from the hypothesized value \(\mu_0 = 95\). We reject the null hypothesis for any mean value outside the “non-rejection interval” \[ \mu_0 \pm F^{-1}_{t(n-1)}(1-\alpha/2) \sqrt{\frac{\sigma^2}{n}} \]

mu_0 <- 95
n <- length(glucose_data)
se <- se_mu_hat(n = n)
margin <- qt(0.975, df = n - 1) * se
upperbound <- mu_0 + margin
lowerbound <- mu_0 - margin

In this case, the non-rejection interval is \([92.959028, 97.040972]\).

Calculate power under a simple alternative

Consider the simple alternative that the true value is actually the estimated mean calculated from the data (i.e. \(98.66\)). Let’s also assume that the known standard deviation is what we estimated from the data.

prob_low <- pt(
  q = (lowerbound - mu_hat) / se,
  df = n - 1,
  lower.tail = TRUE
)

prob_high <- pt(
  q = (upperbound - mu_hat) / se,
  df = n - 1,
  lower.tail = FALSE
)

power <- prob_low + prob_high
print(power)
#> [1] 0.940662

Power as a function of sample size

power <- function(n = 100, null = 95, alt = 98.66) {
  # there's no such thing as fractional sample size:
  n <- floor(n)
  # using the function we wrote earlier:
  se <- se_mu_hat(n = n)
  reject_upper <- ((null + qt(0.975, df = n - 1) * se) - alt) / se
  reject_lower <- ((null - qt(0.975, df = n - 1) * se) - alt) / se
  p_reject_high <-
    pt(
      q = reject_lower,
      df = n - 1
    )
  p_reject_low <-
    pt(
      q = reject_upper,
      df = n - 1,
      lower = FALSE
    )
  p_reject <- p_reject_high + p_reject_low
  return(p_reject)
}
power_plot <-
  ggplot() +
  geom_function(fun = power, n = 100) +
  xlim(c(2, 200)) + # n = 1 is not allowed for t-distribution
  ylim(0, 1) +
  ylab("Power") +
  xlab("n") +
  theme_bw()
print(power_plot)

Simulations

Create simulation framework

Here’s a function that performs a single simulation of a Gaussian modeling analysis:

do_one_sim <- function(
  n = 100,
  mu = mean(glucose_data),
  mu_0 = mean(glucose_data) * 0.9,
  sigma2 = var(glucose_data),
  return_data = FALSE # if this is set to true, we will create a list() 
  # containing both the analytic results and the vector of simulated data
) {
  # generate data
  x <- rnorm(n = 100, mean = mu, sd = sqrt(sigma2))

  # analyze data
  mu_hat <- mean(x)
  sigmahat <- sd(x)
  se_hat <- sigmahat / sqrt(n)
  confint <- mu_hat + c(-1, 1) * se_hat * qt(.975, df = n - 1)
  tstat <- abs(mu_hat - mu_0) / se_hat
  pval <- pt(df = n - 1, q = tstat, lower = FALSE) * 2
  confint_covers <- between(mu, confint[1], confint[2])
  test_rejects <- pval < 0.05

  # if you want spaces, hyphens, or characters in your column names, 
  # use "", '', or ``:
  to_return <- tibble(
    "mu-hat" = mu_hat,
    "sigma-hat" = sigmahat,
    "se_hat" = se_hat,
    "confint_left" = confint[1],
    "confint_right" = confint[2],
    "tstat" = tstat,
    "pval" = pval,
    "confint covers true mu" = confint_covers,
    "test rejects null hypothesis" = test_rejects
  )

  if (return_data) {
    return(
      list(
        data = x,
        results = to_return
      )
    )
  } else {
    return(to_return)
  }
}

Let’s see what this function outputs for us:

do_one_sim()

Looks good!

Now let’s check it against the t.test() function from the stats package:

set.seed(1)
mu <- mean(glucose_data)
mu_0 <- 80
sim_output <- do_one_sim(mu_0 = mu_0, return_data = TRUE)
our_results <-
  sim_output$results |>
  mutate(source = "`do_one_sim()`")

results_t_test <- t.test(sim_output$data, mu = mu_0)

results2 <-
  tibble(
    source = "`stats::t.test()`",
    "mu-hat" = results_t_test$estimate,
    "sigma-hat" = results_t_test$stderr * sqrt(length(sim_output$data)),
    "se_hat" = results_t_test$stderr,
    confint_left = results_t_test$conf.int[1],
    confint_right = results_t_test$conf.int[2],
    tstat = results_t_test$statistic,
    pval = results_t_test$p.value,
    "confint covers true mu" = between(mu, confint_left, confint_right),
    `test rejects null hypothesis` = pval < 0.05
  )

comparison <-
  bind_rows(
    our_results,
    results2
  ) |>
  relocate(
    "source",
    .before = everything()
  )

comparison

Looks like we got it right!

Run 1000 simulations

Here’s a function that calls the previous function n_sims times and summarizes the results:

do_n_sims <- function(
  n_sims = 1000,
  ... # this symbol means "allow additional arguments to be passed on to the 
  # `do_sim_once` function
) {
  sim_results <- NULL # we're going to create a "tibble" of results,
  # row by row (slightly different from the hint on the homework)

  for (i in 1:n_sims) {
    set.seed(i) # sets a different seed for each simulation iteration, 
    # to get a different dataset each time

    current_results <-
      do_one_sim(...) |> # here's where the simulation actually gets run
      mutate(
        sim_number = i
      ) |>
      relocate("sim_number", .before = everything())

    sim_results <-
      sim_results |>
      bind_rows(current_results)
  }

  return(sim_results)
}

sim_results <- do_n_sims(
  n_sims = 1000,
  mu = mean(glucose_data),
  sigma2 = var(glucose_data),
  n = 100 # this is the number of samples per simulated data set
)

sim_results

The simulation results are in! Now we have to analyze them.

Analyze simulation results

To do that, we write another function:

summarize_sim <- function(
    sim_results,
    mu = mean(glucose_data),
    sigma2 = var(glucose_data),
    n = 100) {
  # calculate the true standard error based on the data-generating parameters:
  se_mu_hat <- sqrt(sigma2 / n)

  sim_results |>
    summarize(
      `bias[mu-hat]` = mean(.data$`mu-hat`) - mu,
      `SE(mu-hat)` = sd(.data$`mu-hat`),
      `bias[SE-hat]` = mean(.data$se_hat) - se_mu_hat,
      `SE(SE-hat)` = sd(.data$se_hat),
      coverage = mean(.data$`confint covers true mu`),
      power = mean(.data$`test rejects null hypothesis`)
    )
}

Let’s try it out:

sim_summary <- summarize_sim(
  sim_results,
  mu = mean(glucose_data),
  # this function needs to know the true parameter values in order to assess 
  # bias
  sigma2 = var(glucose_data),
  n = 100
)

sim_summary

From this simulation, we observe that our estimate of \(\mu\), \(\hat\mu\), has minimal bias, and so does our estimate of \(SE(\hat\mu)\), \(\hat{SE}(\hat\mu)\).

The confidence intervals captured the true value even more often than they were supposed to, and the hypothesis test always rejected the null hypothesis.

I wonder what would happen with a different sample size, a different true \(\mu\) value, or a different \(\sigma^2\) value…

Graph the Likelihood

mu_likplot <- 
  ggplot() + 
  geom_function(fun = function(x) lik(mu = x)) + 
  xlim(mean(glucose_data) + c(-1,1) * sd(glucose_data)) +
  ylab("likelihood") + 
  xlab("mu") +
  geom_vline(xintercept = mean(glucose_data), col = "red")

Graph the Log-likelihood

ggplot() + 
  geom_function(fun = function(x) loglik(mu = x)) + 
  xlim(mean(glucose_data) + c(-1,1) * sd(glucose_data)) +
  ylab('log(likelihood)') + 
  xlab("mu") +
  geom_vline(xintercept = mean(glucose_data), col = "red")